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Stability of freely falling granular streams 
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A freely falling stream of weakly cohesive granular particles is modeled and analysed with help of 
event driven simulations and continuum hydrodynamics. The former show a breakup of the stream 
into droplets, whose size is measured as a function of cohesive energy. Extensional flow is an exact 
solution of the one-dimensional Navier- Stokes equation, corresponding to a strain rate, decaying 
like t~^ from its initial value, 70. Expanding around this basic state, we show that the flow is stable 
for short times, 70 1 ^ 1, whereas for long times, 70 1 ^ 1, perturbations of all wavelength grow. 
The growthrate of a given wavelength depends on the instant of time when the fluctuation occurs, 
so that the observable patterns can vary considerably. 

PACS numbers: 83.50.Jf, 47.20.-k, 45.70.-n 



Recent experiments on granular streams have revealed 
many features which are familiar from molecular liquids. 
Somewhat surprising was the observation of clustering 
[TH5 in freely falling dry granular streams which are rem- 
iniscent of the droplet patterns observed in liquids due to 
surface tension. Even though tiny attrative forces could 
be measured and are attributed to van der Waals in- 
teractions or capillary bridges, the observed size of the 
clusters did not agree with the predictions of Rayleigh- 
Plateau. In another set of experiments [H |5| capillary 
waves and their dispersion were measured, allowing to 
deduce a (tiny) surface tension. Exciting perturbations 
of a given frequency and observing their initial growth 
was consistent with the Rayleigh-Plateau analysis. 

In this paper, we model a freely expanding stream of 
weakly cohesive, inelastically colliding grains and simu- 
late it for the parameters deduced from experiment. We 
confirm the observed clustering and determine growth 
rates and drop sizes in dependence on cohesive energy. 
The initial instability is analysed within a continuum de- 
scription, based on the Navier-Stokes equations. Given 
an exact solution of the nonlinear equtions for extensional 
flow, linear stability analysis can be performed and pre- 
dicts nonmonotoneous behaviour as a function of time: 
For short times a finite strainrate stabilises the stream, 
wheras for long times it becomes completely unstable. 

Cohesive forces — We model [6 the grains as hard 
spheres of diameter d. When two particles approach they 
do not interact until they are in contact whereupon they 
are inelastically reflected with a coefficient of restitution 
£. Moving apart, the particles feel an attractive poten- 
tial of range d^. Such an attractive force can be due to 
capillary bridges or van der Waal forces, if the particles 
are deformed in collisions. As the spheres withdraw be- 
yond the distance dcf , a constant amount of energy VFcoh 
is lost provided the normal relative velocity Av of the 
impacting particles is sufficien t to over come the poten- 
tial barrier, /S.v > A'Ucrit = \/2Wcoh/M (where /i is the 
reduced mass), otherwise the particles form a bounded 
state, oscillating back and forth. 



Strainrate — We assume that the particles fall out of 
the container into a vacuum [7] with an intial velocity vq. 
For simplicity, consider a column of n particles leaving 
the hopper sequentially, with a time intervall At = d/vo, 
and ignore collisions for now. We notice that the i^^ 
particle will be accelerated according to Zi{t) = g{t — 
iAt) + '^0- Hence there is an initial velocity gradient, 
which can be computed from 
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In the comoving frame the stream expands freely: the 
particles move with constant velocity, however their dis- 
tance increases. Hence we expect the strain rate to de- 
crease as a function of time according to 
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dv 

dz 1 + 7ot 



(2) 



This will turn out to be important for the stability analy- 
sis: stretching is known [8 to stabilise the flow and hence 
prevent clustering. As we will see below, this is precisely 
what happens for short times, wheras for long times we 
recover the clustering instability, when the strain rate has 
become sufficiently small. 

Simulation — This simple model can be simulated with 
an event driven code [6 , allowing us to consider large 
systems with up to A/" = 10^ particles. We simulate the 
freely falling stream in the rest frame of the stream, im- 
posing a homogeneous velocity gradient 7o = ^ = ^ 
together with a small random velocity in the initial state. 
Given the instantaneous interactions in our simple model, 
the strain rate 70 and the cohesive energy Wcoh are not 
independent parameters: If, e.g., 70 is increased by a fac- 
tor of two and VFcoh by a factor of 4, the particles follow 
exactly the same trajectories, just twice as fast. Hence, 
we only vary the cohesive energy, w := Wcoh/I^coh,exp, 
which is convenienly measured relative to the typical ex- 
perimental value of ref. [2 , i.e. Wcoh,exp = 10~^^ J. The 
cohesive energy, relates to the surface tension F, used 



later, through F ^ Wcoh/d'^ El E]. The remaining pa- 
rameters are chosen, unless specified otherwise, to match 
the typical experimental values, namely the coefficient 
of restitution e = 0.9 , stream's initial volume fraction 
(f) = 0.5, and initial stream radius r^ = 19d. 

Simulation results — Fig. [l] shows snapshots of the 
same system at five different times, demonstrating, how 
the initially straight stream profile develops inhomo- 
geneities which grow in time and finally lead to separate 
clusters. 
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FIG. 1: (color online) Snapshots of the system for different 
times; colors indicate relative motion of adjacent particles; 
frozen areas appear in gray (dark) and yellow (bright) cor- 
responds to areas of ongoing deformations; small droplets of 
size Ndvop < 1000 are ignored for better visibility. See [lOj for 
a movie. 

In the inset of Fig. [2] we plot the mean droplet size 
A^drop as a function of time for two values of Wcoh- After 
a sharp initial decrease due separation of the stream into 
clusters, A^drop reaches a steady state. Its value is shown 
in the main plot for a range of cohesive energies. Scal- 
ing arguments in [3 suggest that the typical length of a 
droplet, rescaled back to its length on the unstretched 
stream, Aq, should scale like the square root of the co- 

1 /2 

hesive energy. Hence, we expect A/drop ex Aq ex W^^^i • 
The solid line is a fit to the data points with an exponent 
f3 = 0.54, confirming the simple scaling arguments. 

The actual shape of the droplet is more difficult to 
capture systematically than its mass, since it continues 
to change slightly even after the droplets have separated. 
Royer et al. [2 characterize droplets by their length Ac 
and width Wc , right before they hit the bottom of the ex- 
perimental setup. They find that droplets' aspect ratios 
Xc/wc always fall in between 1 and 3. Even though the 
droplet formation appears to be surface tension driven, 
these findings preclude the expected Rayleigh-Plateau in- 
stability as a predominant mechanism (which only allows 
aspect ratios > tt). In Fig. [s] we show the simulation re- 
sults for the droplet lengths and width for various Wcoh- 
The most striking feature in this plot is the huge scatter 
in droplet length for a given value of Wcoh- This result is 
at variance with a well defined critical wavelength, cor- 
responding to the fastest growing mode. 



FIG. 2: (color online) Mean droplet size (A^drop) as a func- 
tion of w = Wcoh / Wcoh, exp; data points are results from the 
simulation and the solid line is a power law fit; inset: mean 
droplet size (A^drop) as a function of time for it; = 8 (high final 
value) and w = 1 (low final value); at the grey vertical line 
all systems have reached a steady state, which is used for the 
main plot. 
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FIG. 3: (color online). Length vs. width of individual clusters 



with A^drop > 1000; w ^ 1 (disk), 



3 (triangle), w — 8 



(square), and w — 15 (star); gray lines correspond to aspect 
ratios of 1 and 3; inset: time Tc for an imposed perturbation 
of wavenumber K to grow beyond a given value Ac (1.7, 1.5, 
1.3, 1.1 from top to bottom) 



Previous findings suggest that droplet formation is due 
to an instability, causing small fiuctuations at certain 
wavelengths to grow, while other wavelengths are stable. 
To shed light on this instability, we impose a small undu- 
lation h{z) = To + ecos(Az) in the initial state, follow the 
time evolution of the respective Fourier mode A(A, t) and 
determine the time Tc, it takes the amplitude to grow be- 
yond a certain value Ac, i.e. A(A, rc)/A(A, 0) = Ac. This 
result is shown in the inset of Fig. [3J A fastest growing 
mode can be identified and hardly depends on the choice 
of Ac. In the following section, we study this instability 
in terms of a continuum theory and compare the predic- 
tions to simulation and experiment. 

Continuum theory — To analyse the stability of the 



initially homogeneous stream we use continuum theory 
[8l[TT]. Our starting point are the Navier-Stokes equa- 
tions for the velocity field, iT(r, z;t), in cylindrical coor- 
dinates assuming axial symmetry 



dtv ^ {v - V)v ■■ 
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(3) 



together with the equation of motion for the interface 
r = h{z,t), 



dth^Vzdzh = Vr\r^ 



(4) 



Here p denotes the pressure, p the density and u the shear 
viscosity. These equations have to be solved, subject to 
the boundary conditions, requiring the balance of normal 
and tangential forces at the interface: an = —nVpn 
8it r = h. Here K: is the curvature of the interface, T 
is the surface tension divided by the density and aij = 
—pSij -\- u{diVj + djVi)/p denotes the stress tensor. 

To obtain approximate solutions to the above equa- 
tions, we follow Eggers [8 and assume that variations in 
the radial direction take place on scales small compared 
to variations along the stream. Under these assumtions 
a one dimensional Navier Stokes equation for v = Vz{z^ t) 
has been derived ^ for an incompressible fluid: 
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(6) 



These equations have been studied in various circum- 
stances for molecular fluids [HI [11] . The best known one 
is the Rayleigh Plateau instability, where one expands 
around a state with constant radius and velocity which 
does not apply in the presence of gravity. Jet flow dom- 
inated by viscous effects [121 [E] has also been analysed 
within the above one- dimensional model. Here we con- 
sider instead a freely falling stream [TT in the comoving 
frame. This state is characterized by a time dependent 
velocity gradient, that is constant in space: v{z^t) = 
i^^' Incompressibilty requires h{z^t) = ro(l+7o^)~"^^^- 
These fields solve the above equations exactly^ allowing 
us to do a linear stability analysis by expanding around 
the above solution. 

We introduce a dimensionless position variable Z := 
z/{ro{l + 7o^)) such that Z remains fixed, if z moves 
along with the stream. The Z-dependence can then be 
taken care of by plane waves: ~ exp {ikr^Z). To further 
simplify the notation, we introduce dimensionless time 
r = 7ot and wavenumber K = kr^. We obtain two linear 
equations for h{z^t) — h{z^t) = ex.-p{iKZ)eR{'yot) and 
v{z, t) - v{z, t) = exp {iKZ)joev{jot): 
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(7) 



For clarity of presentation we have set z/ = and hence 
are left with one dimensionless parameter F = F/(rQ7Q). 
The generalisation to finite viscosity is straighforward 
and given in the supplementary material [TF. 

The above eqs. are two ordinary differential equations 
with time- dependent coefficients. This makes the stabil- 
ity analysis complex, because a given wavenumber, which 
is stable at to, can override an initially unstable mode in 
the course of time. Of course the equations can easily be 
integrated numerically. Before discussing the generalised 
eigenvalue problem, we try to extract the qualitative be- 
haviour by inspecting the equations for small and large 
times. For r = 7ot < 1 we expect the initial strainrate 
to have a stabilising effect, in particular for long wave- 
length perturbations. For long times, r = jot > 1, on 
the other hand, the strainrate decays. With the ansatz 
ansatz ey^^H ^ ^^^ one finds for r ^ 1 
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that all wavenumbers are unstable for sufficiently large 
r, because the dominant term is the one involving the 
surface tension (F). To discuss the general case, we com- 



0.2 



Re(Amax) 



0.1 



0.0 ic 




FIG. 4: Real part of the unstable eigenvalue as a function of 
wavenumber and time 

pute the eigenvalues for all K and t, by diagonalising the 
time-dependent matrix of coefficients. The larger eigen- 
value, Amax, - responsible for the instability - is shown 
in Fig. [4] as a function of K and t in the range of values, 
where Amax > 0. The parameters for the initial strain 
rate, 70, and the surface tension, F, are taken from ex- 
periment and the viscosity is set to zero. We observe 
that initially all wavelength are stable for t = 0, the first 
instability sets in at r ?^ 8 and K ^ 15. This wavenum- 
ber is the initially imposed one and has to be scaled 
down by the stretching factor 1 + r, when the stream has 
been stretched up to time r. Furthermore, wavenumbers 
are measured in units of the initial radius of the stream, 
which decreases according to r{r) = ro/(l + r)^/^. Hence 
to obtain the ratio of wavelength to radius at time r. 



when the instabihty occurs, we have to scale wavenum- 
bers according to K{r) = K{l-\- r)~^/^. If one plots the 
eigenvalue versus K{r)^ one observes a ridge at approx- 
imately i^ ^ 1, implying that at each instant of time r 
- while the stream is stretching - unstable modes with 
roughly the wavelength of Rayleigh Plateau are growing. 
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FIG. 5: Example for growth of perturbations, obtained from 
integrating the linearised equations; inset: most unstable 
wavenumber and time of growth to increase by 50% in de- 
pendence on r/Fexp- 

However, we stress that our system is not in a station- 
ary state, but expanding. As a consequence the eigenval- 
ues are time dependent, and whether a perturbation in- 
creases or decreases depends on the time of its occurence. 
In Fig.|5]we show an example of the integrated Eqs. with 
3 modes excited initially {cr = 1 for K = 1,4,10). As 
expected from the eigenvalue analysis, initially all modes 
are stable, then the smallest K starts to grow, but is - at 
later times - overridden by the larger wavenumber s. This 
is meant to illustrate that the observed pattern will de- 
pend on the instant, when a fluctuation with a particular 
wavenumber occurs. 



Variation of the parameters — If the viscosity is in- 
creased, the results remain qualitatively the same, but 
the instabilities occur at later times, because viscosity 
tends to stabilise the flow. If the initial strain rate is 
put to zero, we recover the Rayleigh-Plateau instability. 
Interesting effects are observed by varying the cohesive 
energy. Decreasing the cohesive energy, F < Fexp, the 
initial range of unstable wavenumbers shifts to larger K^ 
i.e. smaller wavelength in agreement with simulations 
(see Fig.|2]). We determine the time Ta{K) which it takes 
an unstable mode of wavenumber K to grow by 50% from 
its initial value for several values of cohesive energy. In 
that way we can deduce the critical wavenumber and the 
time for the stability to occur as a function of F. These 
are shown in the inset of Fig. |5] We clearly observe an 
increase in critical wavelength with F in agreement with 
simulation and experiment. Furthermore for increased F 
the instability occurs earlier. 

Conclusions — We have shown that a stream of gran- 
ular particles falling under gravity is generically unsta- 
ble due to surface tension - even though the Rayleigh- 
Plateau argument does not apply. In the comoving frame 
the stream is freely expanding, implying that the initially 
straight profile is subject to a time-dependent strain rate. 
Linearising the Navier- Stokes equation around this non- 
stationary state, we have shown that the strain rate sta- 
bilises the straight flow profile at short times, whereas 
for long times all wavenumbers are unstable. Since 
we expand around a nonstationary state, the growth 
rate of a given wavelength depends on the time, when 
the corresponding fluctuation occurs spontaneosuly or 
is introduced into the flow. Thus a variety of patterns 
may be observed including behaviour reminiscent of the 
Rayleigh-Plateau instability [5 . 
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